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ABSTRACT 


The  purpose  of  this  work  is  to  present  an  inversion  algorithm  fox 
backscattered  ('"'stacked*'')  seismic  data  which  will  reconstxuct  the  velocity 
profile  in  realistic  earth  conditions.  The  basic  approach  follows  that  of 
the  original  Cohen  and  Blei  stein  paper  (197S>aK^Tn  that  high  frequency 
asymptotics  and  perturbation  methods  are  used.  However,  in  the  original 
paper  the  perturbation  was  rexative  to  a  constant  reference  speed,  whereas 
the  current  work  uses  a  reference  speed  which  may  vary  with  depth.  This 
greatly  enhances  the  validity  of  the  perturbation  assumption  and  hence  the 
inversion  results.  On  the  other  hand,  the  new  algorithm  enjoys  the  same 
economies  and  stability  properties  of  the  original  algorithm,  making  it  very 
competitive  with  current  migration  schemes. 

Four  major  assumptions  are  made:  (|)  the  acoustic  wave  equation  is  an 

A 

adequate  model,  (id:)  stacked  data  has  amplitude  information  worth  preserving 

A 

fairly  accurately,  (ddi)  the  actual  reflectivity  coefficients  can  be 

A 

adequately  modeled  as  perturbations  from  a  continuous  reference  velocity 
which  depends  only  on  the  depth  variable,  and  tiwd  the  subsurface  can  be 

A 

adequately  modeled  as  a  series  of  layers  with  jump  discontinuities  in  the 
velocity  (or  impedance)  at  these  layers. 

\ 

While  the  algorithm  is  particularly  suited  for  data  generated  by  a 
number  of  reflecting  surfaces,  its  validity  for  a  single  reflector  is 
demonstrated  by  applying  the  algorithm  to  Kirchhoff  data  for  a  quite  general 
surface. 


i 


A  key  feature  of  the  approach  of  this  paper  is  the  repeated  application 
of  high  frequency  asymptotic  methods;  both  in  obtaining  the  basic  integral 
equation  describing  the  unknown  velocity  correction,  and  in  the  inversion  of 
this  integral  equation.  Perhaps  a  noteworthy  feature  is  that  the  underlying 
integral  equation  is  in  the  form  of  a  generalized  Fourier  integral  equation; 
and  the  method  for  its  (approximate)  inversion  may  prove  to  be  applicable  to 
a  wide  range  of  such  problems. 


GLOSSARY 


b(r, r’) 

B 

<=1 

c(z) 

D,  E.  F.  G,  H  ( K,  z ) 
f 

g(x,  z;£. w) 
h(x) 

I(r.r') 

K 

k.(K.z) 

n(z)  *  c(0)/c(z), 
r  =  (x,  y,  z) 

R 

Rj 

s.  Sj 
S 

ns<£.-fc>) 

Ds(J,t) 
v(x»  z) 

W 

X  =  (x.  y) 


amplitude  of  Green's  function  (B-l) 
a  basic  amplitude  (14) 
inversion  amplitude  (10) 
velocity  below  reflector  (30) 
reference  velocity  (1) 
various  integrals;  see  Appendix  A 
frequency  (46) 

Green’s  function  for  c(z)  medium  (2) 
cylindrical  surface  (25) 
delta-function-like  integral  (12) 
ray  parameter  (3) 

Jn*(z)-K*  (6) 

the  index  of  refraction  (7) 
cartesian  (2) 

reflection  coefficient  (26) 
reflection  coefficient  (23) 
arclength  variables  (23) 
abbreviation;  see  (29) 
observed  scattered  field  (2) 
observed  scattered  field  (45) 
velocity  (1) 

the  inversion  operator  (10) 
horizontal  cartesian  (1) 
vertical  cartesian  (1) 
unknown  perturbation  in  velocity  (1) 


u(x, z) 
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p 

T  (  K,  Z  ) 

f 

u 


reflectivity  function  (23) 
abbreviation.*  see  (26) 
abbreviation.*  see  (30) 

Dirak  delta  function,*  belov  (12) 
band-limited  distributions  (18) 
jump  in  c  (31) 

cartesians  for  observation  point  at  z  =  0  (2) 

signed  distance  from  r'  (20) 

travel  time  (8) 

phase  function  (14) 

circular  frequency  (2) 


INTRODUCTION 


Carter  and  Frazer  [1984],  and  Blei stein  and  Gray  [1984]  (henceforth  BG), 


present  i 

mversion 

algorithms 

wh  i 

ich  incl 

ude 

the 

effect 

of 

a  stratified 

reference 

velocity, 

c(z)  .  Those 

papers 

did 

not 

address 

the 

question  of 

obtaining 

accurate 

values  of 

the 

reflecti 

on  * 

coef f i 

icient,* 

this 

is  the  issue 

treated  here.  Thus,  in  the  language  of  Bleistein,  Cohen  and  Hagin  [1984], 
(henceforth  BCH) ,  the  earlier  algorithms  provided  structural  inversions,  the 
location  of  the  sub-surface  layers,*  whereas  the  present  algorithm  also 
provides  an  accurate  estimate  of  the  reflectivity  function,  which  depicts 
the  reflectors  and  provides  an  estimate  of  the  reflection  strengths  across 
the  layers. 

Since  we  employ  a  perturbation  assumption  (the  "Born  Approximation"),  the 
constant  reference  speed  inversion  first  described  in  Bleistein  and  Cohen 
[1979a]  and  reviewed  in  BCH,  is  often  not  adequate  at  depth.  Although 
recursive  use  of  the  algorithm  is  possible  and  although  the  results  can  be 
significantly  enhanced  by  suitable  pre~  or  post-processing  (e.g.,  see  Hagin 
and  Cohen  [1984]),  extension  of  the  perturbation  method  to  a  stratified 
reference  profile  is  highly  significant.  It  is  far  more  likely  that  the 
actual  velocity  function  can  be  well  approximated  by  a  stratified  reference 
velocity  than  by  a  constant  one,  which  in  turn  enhances  the  validity  of  the 
perturbation  assumption  and  the  inversion  results.  See  BG  for  further 
discussion  of  this  point. 

The  algorithm  presented  here  has  the  same  structure  as  the  BG  algorithm 
and  hence  it  can  be  expected  to  exhibit  the  same  stability  and  economy.  In 
particular,  we  note  that  the  processing  times  for  this  algorithm  with  depth- 


dependent  background  velocity  will  be  comparable  to  those  for  a  constant 
background  k-f  algorithm.  In  addition,  we  shall  show  below  that  the 
algorithm  can  be  expected  to  be  quite  robust  even  when  the  "small" 
perturbation  assumption  is  violated. 

A  key  feature  of  our  approach  to  this  problem  is  repeated  application  of 
high  frequency  asymptotic  methods  to  obtain  an  inversion  formula  valid  in 
the  high  frequency  regime.  Discussion  of  the  motivation  and  justification 
for  such  high  frequency  approximation  may  be  found  in  BG  and  BCH.  In 
particular,  we  shall  use  a  ray  theoretic  Green's  function  in  formulating  our 
basic  integral  equation,*  equation  (2)  below.  A  similar  approach  was 
presented  in  Clayton  and  Stolt  [1981],  Moreover,  since  the  resulting 
integral  equation  cannot,  in  general,  be  inverted  exactly,  we  also  use 
asymptotic  methods  for  this  task.  This  is  carried  out  in  the  next  section. 
It  is  perhaps  noteworthy  that  this  integral  equation  can  be  viewed  as  a 
generalized  Fourier  integral  equation;  and  hence  the  method  of  inversion  may 
prove  to  be  of  interest  outside  the  present  context.  In  the  subsequent 
section  the  resulting  algorithm  is  verified  (asymptotically)  for  a  quite 
general  class  of  reflecting  surfaces.  A  short  section  on  computational 
considerations  is  included.  Finally,  many  of  the  detailed  calculations  are 
carried  out  in  Appendices  A-E. 


HIGH  FREQUENCY  INVERSION 


Here  we  describe  the  formalism  for  determination  of  an  asymptotic 
inversion  operator.  We  employ  the  same  wave  equation  model  as  described  in 
detail  in  BCH.  If  v  is  the  velocity  in  the  wave  equation,  we  set 


v  (x, z) 


c*  (z) 


1  +  a(  x 


■ 

,  z) 


( x,  y) 


(1) 


We  will  also  use  r  =  (x.z)  =  (x,y,z).  Here  c(z)  is  the  known  stratified 
reference  velocity,  while  a(r)  is  the  desired  perturbation  correction  to  the 
actual  velocity.  Furthermore,  we  retain  the  assumption  of  backscatter 
("stacked”)  data.  In  this  case  the  basic  integral  equation  for  a(r)  is  (cf. 
BCH,  equation  (8)): 


i(  r) 


UgfJ.-w)  =  «  |  |  |  d*r  — - - gl(r;£,*fc>) 

c  (z) 


£  =  U»n> 


(2) 


where  all  unmarked  integral  signs  are  over  (-“>,<=).  Here  u  denotes  the 
backscattered  field  at  the  location  £  =  (£»n)  on  the  observation  plane,  z  = 
0,  and  g  (the  "incident  field”)  denotes  the  Green's  function  corresponding 
to  the  stratification,  c(z)  .  In  contrast  to  the  constant  background  case,  g 
cannot  be  determined  exactly/  we  must  use  the  high  frequency  assumption. 
Fortunately,  this  assumption  is  completely  justified  on  the  geophysical 
exploration  scale  and  has  long  been  used  to  simplify  processing  formulas 
even  when  it  was  possible  to  derive  wide  band  analytic  results  (see  BCH). 
We  use  J.  B.  Keller's  [1978]  ray  method  formalism  (see  also  Bleistein 
[1984]),  which  is  the  multi-dimensional  analogue  of  the  WKB  method  to  obtain 
a  parametric  representation  of  g  (see  Appendix  B) : 


-  3 


(3) 


!*  1  ■■  1  *  V  L" 


i(ur(  K,  z) 

g  ~  -  ■■  6  - - - -  . 

4n  ^  kj(K,  0)  k3(K,  z  )  E(K,  z  )  H(K,  z) 

Here,  if  we  introduce  the  transverse  distance, 

I*  "  &|  «  (4) 

then  the  parameter,  K,  in  (3)  is  defined  as  a  function  of  |x  -  £|  and  z  by 
the  ray  egnation: 

U  -  £|  =  KE(K,z)  .  (5) 

Further,  the  quantity  kj(K#z)  is  given  by 

k(K.z)  =  Jn*(z)  -  K*  ,  (K*  <  n*(z))  (6) 


where  in  turn,  n,  the  index  of  refraction,  is 


n(z) 


_  c  ( 0) 
c(  z)  ' 


(7) 


The  travel-time,  r,  is 


T 


=  'cToT  G(K>  z) 


nNpdt 
k  (K,  f) 

0 


(8) 


and  finally,  the  quantities  E  and  H  are  likewise  integrals  involving  n  and 
k^ .  These  integrals,  as  well  as  others  that  occur  subsequently,  are  defined 
in  Appendix  A.  There,  we  also  derive  some  needed  relations  involving  these 
quantities . 


-  4  - 


We  shall  not  need  the  extension  of  k  in  (6)  to  the  range  n  (z)  <  K 

i 

because  our  Fourier  transforms  are  integrals  over  real  wave  numbers  only. 


Thus  our  task  is  to  invert  the  integral  equation. 


Ug(£»  w) 


,  .  2  i  o*t  (  K,  z ) 

t  a(r)  e 

d  r  2,  :  k  ( K,  0)  k  (K,  z  )  eTK,  Z  )  H(K,  z) 

c  ( z )  1  3 


(9) 


for  a(r)  in  terms  of  the  data,  u  (£, w).  Again,  the  ray  parameter 
K  is  de f i ne d  by  ( 5 ) . 


Since  the  phase  in  (9)  resembles  that  of  a  Fourier  transform,  we  are 
motivated  to  seek  an  asymptotic  inversion  operator  W  of  the  form: 

Wtu^.o))]  (r')  -  {  J  d*$  J  dw  B(r',£)  ,  (10) 

where  the  amplitude  B(r',£)  and  value  K’  must  be  determined.  Condition  (5) 
above  suggests  that  K'  satisfies  Jx'  -  £|  =  K'E(K',z).  In  (10)  we  have 


introduced  primes  to  avoid  confusion  with  the  integration  variables  in  (9). 
Applying  W  to  both  sides  of  (9)  and  writing  out  the  right  side  explicitly  we 


1 


a(r') 


WfugHr') 


16n* 


d*  r 


a  (r)  4  exp{2iw[-c(K,  z)-t(K',  z')]} 

— -  a)  B{  r,  £  ' ) - 

c  (z)  kI(K,0)kj(K,z)E(K.z)H(K.z) 

(11) 

:  |  |  |  d  r  a(r)  I(r,r')/  where 


I(i,r’) 


v  2i0)[r(K,  z)-t(K',  z  ')] 
B(r ',  £/  e 

c’  (x)k,(K.0)k,(K.  z)E(K.z)H(K,z) 


(12) 


Clearly,  if  (11)  is  to  hold,  then  I(r,r')  must  represent  the  3-dimensional 
delta  function,  8(r-r*).  Hence  our  task  is  to  find  B  =  B(r',£)  and  K',  so 
that  this  is  the  case. 


First  a  comment  about  the  assumed  form  of  B,  i.e.,  B(r',£)  where 
r  =  ( x, y,  z)  and  £  =  (£,ii).  In  attempting  to  invert  (9)  for  a(r’)  one  would 
generally  expect  B  to  also  involve  w.  However,  our  experience  with 
canonical  problems  (e.g.  c(z)  =  constant)  suggests  that  B  is  independent  of 
to,  and  the  work  to  follow  confirms  this. 


In  (11)  and  (12)  we  must  acknowledge  that  in  Geophysical  applications  <d 
is  confined  to  the  high  frequency  regime  (e.g.  see  BCH  for  details).  Hence 
we  may  evaluate  the  asymptotic  (large  <o)  approximation  to  the  £  integral  in 
(12)  by  stationary  phase.  For  convenience  we  will  think  of  w  as  the  "large 
parameter”  in  our  asymptotics  to  follow.  More  precisely,  it  is  important 


that  the  dimensionless  parameter  2uiL/c(0)  be  "large*,  where  L  is  a  typical 


length  scale  of  the  problem.  This  is  discussed  more  fully  in  BG.  The 
details  of  the  stationary  phase  calculation  are  carried  out  in  Appendix  D. 
The  conclusions  are: 

1)  K'  =  K.  Geometrically,  this  stationarity  condition  says  that  for  given  r 

and  r',  £  =  £  is  chosen  so  that  r,  r'  and  (£  ,0)  lie  on  the  same  ray. 
s  s 

See  Fig.  1.  Since  K'  =  K  we  will  use  only  K  hereafter.  Hence,  along 
with  (5),  we  now  have 

|x*  -  £|  =  K  E(K,  z  *)  .  (13) 


2)  The  asymptotic  approximation  to  (12)  is  as  follows  where  E'  =  E(K,  z'),  H' 

=  H(K,  z')  and  £  symbolizes  that  the  stationarity  condition,  K’  =  K,  is 
s 

to  be  applied: 


Hr,  r ' ) 


1 

16  TT* 


I 


dm  to 


b(r,  r') 


2iti>[T(K,  z) 
e 


r(K,  z')]j 


where 


w  _  B  nc(0)  sgn(z'-z)  (E'H'/EH)/ 

b(r.r')  - S - 

c  (z)  k  (K.0)  k  (K,  z)  [(E  -  E')(H  -  H’)]  ' 


(14) 


This  expression  for  I(r,r’)  will  simplify  dramatically  after  we  perform 


the  to  integration.  At  this  point  we  have 


I(r,r') 


b(r.r')  f 

16  n*  J 


da)  ib>  e 


2iu>[x(K.z)  -  x(K.z')] 


If  we  denote  x'  =  x(K.z’),  those  familar  with  distributions  will  recognize 
the  above  w  integral  as  the  distribution 


f  s‘(x-x')  ^  8(x-x  * ) 


where  5  is  the  standard  delta  function.  The  distribution  in  (16)  has  the 
sifting  property  that,  for  any  differentiable  f(x). 


f(x)  6'(x-x  )dx  =  -  f'(x  ) 
o  o 


We  will  use  this  property  shortly.  We  now  have 


I(r.r’)  ~ 


b(r,  r')  , 


8b(x-x') 


where  the  subscript  B  serves  as  a  reminder  thst,  in  application  of  these 
results  the  frequency  u  will  actually  be  band-limited.  Hence  the  w  integral 
in  (IS)  is  proper  and,  consequently,  the  resulting  distribution  6'  is  in 
fact  a  band-limited  version  (i.e.,  a  regular  function).  These  matters  have 
been  discussed  in  Cohen  and  Bleistein  [1979b],  Bleistein  [1984]  and  BCH 

f 

[1984].  We  will  use  the  notation  and  for  this  reason,  but  proceed  as 


if  all  the  action  in  I(r,r')  takes  place  as  r  -)  r*. 

If  one  tarns  to  b(r.r')  defined  by  (14)  and  refers  to  the  definitions  in 
Appendix  A,  it  is  easy  to  show  that  as  r  approaches  r’> 


B  kj (K. x  * ) 

~  32  c(z')kj(K/6)(z-'-zT 


Using  this  in  (18)  gives 


B  k  (K,  x ') 

Kr.r') - - -  ' ) 

32c( z ' )  k} (K,  0)  (z'-z) 


(19) 


Notice  the  singularity  (z'-z)  ,  as  z  ->  z' ,  in  this  form  of  I.  This  is 

precisely  the  amount  of  singularity  needed,  when  combined  with  6'(t-t’),  to 
produce  the  required  3-dimensional  distribution  8( r-r ' ) •  Recall  that  as 
x  ->*'*  £  ->r'  along  the  common  ray  and  x  =  t(K,  z)  t(K,  z’)  =  x'.  In 


Appendix  E  the  following  is  established/  if 


S(t“t ')  —  c  ( z  * )  6  ( p )  =  c  ( z  * )  6  '  ( p) 


n(z') 


z-»z'  z-z'  kj  ( K,  z  ' ) 


We  use  these  two  results  to  rewrite  the  singular  portion  of  (19)/  dropping 


the  subscript  B  for  now. 


6  '  (t-t ' ) 


n(  z  ')  6  ’  ( p) 


We  now  verify  that  -8*(p)/p  behaves  like  2n  8(r-r’)  by  integrating  it 


against  an  arbitrary  differentiable  function  f(r)  over  an  e-sphere  centered 


at  r'. 


Ill 


*tt  e 


dV  f (r) 


S  '  ( p) 


=  j"  d0  sin0  |  d0  |  dp  p*  f(i 


8  *  ( p) 


|r-r ’ |-e 


0  0 


an  e 


dO  si 


n  0  d*S  dp  [p  f(r)l 


0  -e 


=  2n  [p  f(r)l 


2n  f(r') 


.  k  .  .  .  *  '  -  *  -  -  -  "  *  -  !•  ”  *  *  ~  "  *  -m  -  ^  I  W-  »•-  .  W 


Note,  in  the  second  step,  the  compensating  change  of  limits  in  the  0  and  p 
integrals.  In  the  last  step  we  used  the  sifting  property  (17)  of  6*(p).  In 
the  sense  of  3-dimensional  distributions,  we  have  established 


6 '  ( p) 
-P 


2n6(r-r') 


Using  this  and  (20)  in  (19),  and  returning  to  our  6^  notation,  we  have  at 
last 


B  n  c(0)  ,  ,  ,, 

I(*#i  )  ~  16k  (K,0)  8B(i  - 


For  this  to  represent  6_ ( r-r 1 )  we  are  led  to  define  amplitude  B  as 

B 


B  = 


16  k( (K, 0) 
— — 


16  Jl-K* 

n^TCT 


(21) 


This  completes  the  definition  of  our  inversion  operator  V  in  (10). 
Applying  W  to  (9)  gives  our  inversion  for  a(r'). 


-2ia>r(K,  z ') 

e 


Ujj  (J,  w) 


-<£•’  -  1 1  *’«  f 


(22) 


As  discussed  in  BCH,  once  we  surrender  knowledge  of  the  low  frequency 
input  information,  we  cannot  obtain  output  trend  information.  It  is  to  be 
hoped  that  (by  iteration  if  necessary)  our  c(z)  reference  velocity  is  an 
adequate  approximation  of  the  trend  to  the  depths  of  interest.  What  we  can 
obtain  from  band-limited  information  is  a  perturbation  correction  consistent 
with  the  model  of  jumps  across  a  series  of  interfaces.  We  determine  the 
approximate  location  of  these  interfaces  as  well  as  the  approximate  value  of 
the  reflection  coefficient  at  the  interfaces.  This  information  is  summed  up 
in  the  reflectivity  function. 


P  =  l  R.8„( s . )  (23) 

J  °  J 

th 

where  s  is  a  (local)  arclength  variable  measured  normally  from  the  j 

j 

interface  and  R  is  the  normal  reflection  coefficient  of  that  interface. 

j 

Clearly,  knowledge  of  p  is  equivalent  to  knowledge  of  reflector  location  and 
the  normal  reflection  coefficient  (see  equation  (44)  below).  In  turn  the 
latter  allows  direct  computation  of  the  jump  in  c  across  the  reflector. 


According  to  the  theory  developed  in  Cohen  and  Bleistein  (1979b)  and 


reviewed  in  BCH,  we  can  obtain  p  from  a  by  inserting  a  factor  of  i»/2c(z)  in 
(22)  to  obtain  (after  dropping  the  primes): 


Verification  of  Algorithm 


In  order  to  verify  the  feasibility  of  the  algorithm  in  (24),  we  first 
obtain  an  expression  for  the  Kirchhoff  representation  of  data  •  Such  data 
employs  the  high  frequency  assumption  (by  using  the  multi-dimensional  WKB 
representation  of  the  incident,  reflected  and  transmitted  fields),  but  does 
not  make  the  Born  approximation  of  small  reflection  coefficient.  Thus,  if 
we  can  show  asymptotic  equality  in  (24)  for  such  data,  our  algorithm  is 
likely  to  be  quite  robust  for  large  contrast  interfaces. 


It  remains  to  decide  on  the  surfaces  to  use  in  computing  the  Kirchhoff 
approximation  to  u  .  For  ease  of  presentation,  we  carry  out  our 

iJ 

calculations  for  the  general  cylindrical  (i.e.  y  independent)  surface: 


z  =  h(x) 


With  a  little  more  effort  the  verification  can  be  done  for  a  quite  general 
reflecting  surface. 


In  Appendix  C,  we  show  that  the  back  scatter  at  £  from  (25)  has  the 
Kirchhoff  (high  frequency)  representation: 


u^(£,<d)  ~  2iw 


d2X  J 1 +h ' 1 ( X )  y  R  S  e 


2iurt(  K»  X) 


subject  to 


ll  "  l\  «  KE  (K,Z)  »  X  s  h ( x)  . 


Here  we  have  used  bars  to  distinguish  the  spatial  variable  from  the  output 


variables  in  the  inversion  formula  (24)  and  have  also  introduced  the 


quantities : 


[  IIZILL!-  -  k  (i.a)  1 
L  EU,  z)  *  J 


c(o)  n+h' 


16n  ka(K.O)k,(K.Z)E(K.Z)H(K.Z) 


and  the  (non-normal)  reflection  coefficient. 


Y_Yi  s  1  1 

R  =  •  rt  =  sgn(y)  Y  +  T  "  — 

i  ,  Cl  c 


In  turn,  ct(z)  denotes  the  actual  velocity  below  the  reflector,  I  =  h(f). 


that  is. 


c  (z)  =  c(z)  +  Ac 
x 


where  Ac  is  not  accounted  for  by  the  stratified  reference  profile  c(z). 
Obviously  determining  R  is  tantamount  to  determining  Ac  or  Cj(*) • 
Putting  data,  Ug  given  by  (26),  into  (24)  leads  to 

P(^  ~  ncHrcTz)  I  du  0)1  II  d*S  II  d*5 


[77  Jl+h'*  yRS  e2i«rT(K.l)  -  t(K,z)l 


Our  goal  is  to  show  that  the  output  of  this  expression  does  indeed  represent 
the  reflecting  surface. 


5 


We  now  carry  out  stationary  phase  in  x,  y,  £,  r\  (see  Appendix  D  for 


details)  and  obtain  the  stationarity  conditions: 


y  =  tj  =  y 


(33) 


x  -  5  =  -sgn(h')  K  E(K,z) 


(34) 


and 


x  -  %  =  -  sgn(h')  K  E(R,  z) 


(35) 


These  imply: 


K  =  f  =  n(z>Jh>l  , 

l|l+h** 


(36) 


kjiK.!) 


n(z) 


(37) 


where  again  z  =  h(x) . 


Geometrically,  these  conditions  confirm  the  cylindrical  nature  of  our 
reflector  and  show  that  the  output  point  (x,  z)  lies  on  a  specular  ray. 
Furthermore  (33-35)  yield 


Y  - 


(38) 


which,  in  turn,  imply  that  R  reduces  to  the  normal  reflection  coefficient. 


Completing  the  stationary  phase  analysis  (again  see  Appendix  D)  we  find 


that 


p(r)  ~ 


16nn2(z)  l|l+h'*  R  S 
1  n 


c(0)  J  |det (# . . ) 

ij 


I 


dot  e 


2ioi[T(K,z)  -  x(K,z)] 


(40) 


Here, 


1 

(1+h'V 

n(7)h" 

r  i  _ 

l 

E(K,  z)E(K,  z) 

H(K,  z)H(f, z) 

l|l  +  h'* 

L  H(i.  z) 

H(K,z) 

(41) 


However,  the  final  integration  in  u»,  yields  a  delta  function  whose  argument 
can  be  transformed  (using  the  relationships  in  Appendix  A)  to  arclength 
along  rays: 


J  dw  e2ib,t-c(t*I)  "  T(t'Z)l  =  nSgtr (K,  z)  -  x(K,z)l 

=  itc(0)  6BCG(K,  7)  -  G(K,  z) ] 
kj(K,  z) 

=  7TC  ( 0)  -  6D[  z  -  z3 

n* ( z)  B 

nc(O)  •  |  t  6p[s(z)-s(z) ] 

n(z)  B 


(42) 


Here  the  last  equality  involving  the  ray  arclength  variable  follows  from 
(A-8)  and  Appendix  B.  We  may  note  that  since  z  =  z  when  the  delta  function 
"acts,"  the  stationarity  conditions  (33-35)  imply  that  the  output  point 
(x,z)  coincides  with  the  specular  point  (x,Z)  =  (1,  h(X)).  Furthermore, 


when  z  =  T  the  second  term  in  the  large  square  brackets  in  (41)  drops  out 
and  hence  using  the  definition  (29),  w*  *  e  that 


J|det(#..)  16  n  (1+h'  )k,  (K.O)^  (K,z) 


Hence  (40)  reduces  to 


p(r)  ~ 


R  n(Z) 


|l+h '  kj  (I,  z) 


6[s(z)-s(z)] 


*  R  6[s(z)-s(2)l 
n 


Here,  the  last  step  follows  from  (37). 


In  summary,  when  the  Kirchhoff  data  u§  for  a  single  reflecting  surface 
Z  =  h(x)  is  put  into  (24)  and  the  computations  are  carried  out 
asymptotically,  the  inversion  algorithm  faithfully  reproduces  the  surface. 
Moreover,  if  a  linear  combination  of  such  data,  representing  several  such 
surfaces,  were  inserted  into  (24)  then,  in  principle,  the  algorithm  could 
reproduce  an  appropriate  sum  of  responses  as  in  (23).  However,  since  the 
background  c(z)  would  presumably  not  be  exact  beneath  the  first  reflector, 
there  would  be  some  distortion  introduced  into  the  second  and  subsequent 
reflectors.  This  issue  was  discussed  in  Hagin  and  Cohen,  1984. 


REMARKS  ON  DATA  PROCESSING 


The  algorithm  for  the  reflectivity  function  derived  in  the  previous 
sections  is 


(J(x.z)  ~ 


8i 


nc(0)  c(z) 


* 


-2 i ci»r ( K,  z) 


d  £  di»  (d  e 


(dt  Us(t,^)  e1 

A 


U)t 


(45) 


|x  -  £|  =  KE(K,  z) 


Here  Ug  is  the  backscatter  time  data  observed  on  an  areal  array.  For 
actual  data  processing,  it  is  convenient  to  "fold*  the  unphysical  negative 
frequencies  onto  the  positive  ones  by  replacing  <»  by  -to  on  the  interval  (- 
<*>,0).  At  the  same  time,  we  introduce  the  physical  frequency  variable 
(measured  in  Hz) : 


f 


(46) 


and  explicitly  acknowledge  the  bandlimiting  by  introducing  F(f),  a  tapered 
high  pass  filter.  After  these  changes,  we  have: 


•  Im  |  df  f  F(f)  e  4nifr(K'z) 
J0 


(47) 


•  f  dt  U_(t,£)  e2rift 
J0 

|x  -  Si  -  KE  (k.z) 


la  practice,  areal  observations  are  often  not  available  and  instead  only 
a  linear  set  of  data  is  used.  In  this  case  we  cannot  hope  to  reconstruct  a 
three  dimensional  image  of  the  subsurface  and  instead  seek  a  two  dimensional 
slice,  p(x,0, z)  ■  p(x, z),  consistent  with  the  data  available.  Since  the 
data  is  now  independent  of  T), 

0s(t.£)  =0s(t.5.0)  «Us(t,£)  ,  (48) 

and  we  may  carry  out  an  additional  stationary  phase  calculation  in  tj.  The 
stationarity  condition  is 


n  =  y 


(49) 


and  the  analogue  of  (45)  is  found  to  be  (see  Appendix  D  for  details): 


p(x,z)  ~ 


fncTo)  c(z) 


j  dl;  ( 1-K  )E(K, 


z) 


I,  r~. —  — 2ib>t(K,  z) 

dw  lib)  e 

00 

dt  Us(t,5)  e11 

*  A 


jx-?|  =  K  E( K, z)  , 


while  (47)  becomes: 


P(x,z)  ~ 


c(0)  c(z) 


d£  i|  ( 1-K  )  E(K,  z) 


00 

(Re  -  Im)  J  df  F( f ) 


-4nifr(K, z) 


j"  dt  Us(t,^)  e2nift  , 


|x-$l  =  K  E (K,  z)  . 


The  basic  concepts  of  reducing  (SI)  to  a  computer  code  are  the  same  as 
those  discussed  in  BG  for  the  algorithm  presented  there.  Briefly,  the  t  and 
f  integrals  are  performed  routinely  using  an  efficient  FFT  algorithm.  The 
main  complication  in  (51)  lies  in  the  expressions  E(K,z)  and  x(K,z),  both 
being  integrals  defined  by  (A-4)  and  (8)  respectively.  This  is  a  bit  subtle 
in  that  the  parameter  K  (see  Appendix  B)  can  be  viewed  as  determining  the 
starting  angle  for  a  ray  connecting  the  surface  point  (€,0)  to  data  point 
(x,  z)  in  (51).  Therefore,  for  a  given  offset  |x-£|,  K  is  defined  by  the 
implicit  relation  |x-?|  =  KE(K,  z) .  In  computation  this  issue  is  handled 
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quite  efficiently  by  two  tables  for  evaluating  t(K,  z)  and  the  amplitude 
(involving  E  in  (51))  as  functions  of  |x-£|  and  z. 

The  computation  time  of  the  resulting  algorithm,  as  pointed  out  in  BG, 
is  comparable  to  a  standard  k-f  migration  algorithm  with  constant  reference 


CONCLUSIONS 


We  have  presented  the  derivation  of  an  inversion  algorithm  for 
backscattered  ("stacked")  seismic  data.  We  made  four  major  assumptions: 
(i)  the  acoustic  wave  equation  is  an  adequate  model,  (ii)  stacked  data  has 
amplitude  information  worth  preserving  fairly  accurately,  (iii)  the  actual 
reflectivity  coefficients  can  be  adequately  modeled  as  perturbations  from  a 
continuous  reference  velocity  which  depends  only  on  the  depth  variable,  (iv) 
the  subsurface  can  be  adequately  modeled  as  a  series  of  layers  with  jump 
discontinuities  in  the  velocity  (or  impedance)  at  these  layers. 

The  last  assumption  is  unavoidable  given  the  nature  of  the  high  pass  (on 
the  exploration  scale)  data  collected  in  the  field.  The  third  assumption  is 
inherent  in  our  approach  although,  as  pointed  out  above,  the  algorithm  can 
be  expected  to  be  robust  even  when  this  assumption  is  violated.  Also  the 
algorithm  presented  here  represents  a  considerable  improvement  over  earlier 
algorithms,  such  as  Cohen  and  Bleistein  [1979a],  which  perturbed  from  a 
constant  reference  velocity. 

On  the  other  hand,  weakening  of  the  first  two  assumptions  seems  eminently 
feasible  and  we  hope  to  apply  the  techniques  expanded  in  this  article  to 
both  inversion  of  offset  data  ("inversion  before  stack")  and  to  equations 
which  more  accurately  describe  the  wave  propagation  in  the  earth. 


It  is  already  known  (see  BG)  that  algorithms  with  the  structure  of  the 
one  presented  here  are  numerically  stable  and  are  computationally  efficient 


relative  to  other  seismic  data  processing  algorithms. 
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APPENDIX  A 


NOTATIONS  AND  IDENTITIES 


We  define 


n(z) 


_  c  ( 0) 
c(  z) 


kf  (K,  z) 


and  the  integrals. 


D(K,z) 


E(K.z) 


F(K,z) 


G(K,  z) 


c(0)r(K, z) 


H(K,z)  = 


n  (Hi 


k*(K,t) 


( A-7) 


s(K, z)  = 


n(  0  d* 

17 


( A-8) 


Similar  quantities  occur  in  the  r-p  theory,  see  Diebold  and  Stoffa 
[1981] .  Among  the  many  relations  which  link  these  quantities,  we  cite  below 
those  that  are  useful  in  carrying  out  the  calculations  presented  in  this 
paper  and  its  appendices. 


First  of  all,  from  (A-2)  it  follows  that 


D  +  K  E  =  G  , 


(A-9) 


E  +  K  F  =  H  . 


(A-10) 


Next  we  cite  the  k  and  z  partial  derivatives  of  D,  E  and  G  which  follow 
respectively  from  use  of 


(A-ll) 


and  the  Fundamental  Theorem  of  calculus 


D_  =  -  KE  ,  D  =  k  (K,z)  , 


(A-12) 


APPENDIX  B 


THE  STRATIFIED  MEDIA  GREEN'S  FUNCTION 

Using  Keller's  [1978]  "ray  method",  developed  in  the  1950's,  we  seek  a 
high  frequency  approximation, 

g(x, z)  ~  A(x, z)  e  -  ,  x  =  { x, y)  (B-l) 

which  asymptotically  satisfies  the  Helmholtz  equation, 

a 

V*g  +  — ^ —  g  =  -6(x-£)  8(z)  ,  i  =  U.q)  .  (B-2) 

c  (z) 

To  complete  the  specification  of  g,  we  insist  that  it  behave  like  the  free 
space  (i.e.  constant  c)  Green's  function  as  the  field  point,  (x, z), 
approaches  the  source  point,  (£,  0).  This  entails  the  conditions, 

t  ->R/c(0)  .  A  (B-3) 


as  R  ->0,  where 


R2  =  |x-£|*  +  z*  .  (B-4) 

We  substitute  (B-l)  into  (B-2)  and  separately  equate  the  coefficients  of 
w*  and  (o  to  zero  (this  is  the  high  frequency  approximation)  giving  rise  to 
the  eikonal  equation, 

k  •  k  =  n*(z)  ,  k  =  c(0)  Vt  ,  n  =  (B-5) 

and  the  transport  equation. 
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2  k  •  VA  +  (V  •  k)  A  =  0 


(B-6) 


The  former  equation  can  be  solved  by  the  method  of  characteristics  (see 
Bleistein,  1984)  which  reduces  the  problem  to  the  solution  of  a  system  of 
ordinary  differential  equations.  The  first  of  these  equations  are: 

^r=?  •  TS-  K  :  Is  <kx*V  '  k  =  (B-7) 

dK  dk 

-f-  =  0  ,  -= - -  n  n'  ( z)  (B-8) 

do  do 

which  define  the  rays i  o  being  the  ray  parameter.  The  source  term  in  (B-2) 
makes 


x( 0)  =  £  .  z(0)  =  0  .  (B-9) 

a  natural  choice  as  initial  data  for  (B-7)  .  The  data  for  k(0)  consists  of 
an  arbitrary  unit  vector  (cf.  (B-5),  noting  that  n(0)  =1).  To  be  specific 
we  choose 


K(0)  arbitrary  ,  k^O)  =  )|l-K 
where  we  have  introduced 


(B-10) 


K  =  |f|  =  j  k*  +  k*  .  (B-ll) 

From  (B-10)  and  the  fact  that  we  are  not  considering  turned  rays  here,  it 
follows  that  K  =  (k1#  k2)  can  be  viewed  as  the  direction  numbers  of  the  rays 
initiating  from  (J;,  q,  0). 


30  - 


Using  (B-6)  and  (B-7)  we  find 


dr  _  n 
da  c ( 0) 


2  ^  *  (V  •  k)  A  =  0 


(B-12) 


(B-13 ) 


which  together  with  the  data  (B-3)  completes  the  specification  of  g, 
expressed  in  (B-l),  in  terms  of  a  system  of  ordinary  differential  equations. 


Proceeding  to  the  analysis,  we  first  note  that  by  (B-8),  K  =  K(0),  so 
henceforth  we  simply  write  K  for  this  constant  vector.  Then  the  eikonal 
equation  ( B— 5 )  gives  us 


k3(K,z)  =  )|n*  -  K*  . 


(B-14) 


Then  (B-7)  yields 


(B-15) 


«*  -  K* 


x  -  l  =  K  E(K,  z) 


where  E  is  defined  by  (A-4) .  Similarly,  we  find 


T  "  THJT  G(K'z) 


(B-16) 


(B-17) 


where  G  is  defined  by  (A-6). 


If  we  introduce  the  ray  J acobian, 

9 (xj  z) 

J  =  d(K,0) 


and  use  the  fact  that 


V 


k  . 


then  the  transport  equation  (B-€)  can  be  recast  as 


A 


constant 


Thus 


A 


lia  A  Jj 
R-4  0 


which  by  (B-3)  implies  that 


A  =  ——  lim 

J7R-.0 


JL 

4nR 


(B-l 8) 


(B-19) 


(B-20) 


(B-21) 


(B-22) 


We  now  indicate  how  to  obtain  the  partial  derivatives  which  are  the 
elements  of  J. 


First  we  take  the  reciprocal  of  the  second  equation  in  (B-7)  and 


integrate  to  get 


j  k}(K,  z')  “ 


We  think  of  z  =  z(K, a),  differentiate  this  expression  for  o  and  use  do/dkj  = 
0  to  obtain 


9z 

3k. 


-  k.k  F 
i  1 


(B-23) 


Similarly  from  (B-16),  and  using  (B-23),  one  obtains 


=  E  8  .  j  =  1.2 

J 


(B-24) 


Since  3x/3 a  and  dz/da  are  given  directly  by  (B-7),  we  are  now  able  to  form 
J.  A  short  calculation,  involving  the  use  of  (A-10)  yields 


J  =  k  EH 

i 


(B-25) 


It  is  easy  to  show  that  as  R  ->  0  (equivalently  a  ->0). 


TV 


(B-26) 


la  the  seme  limit,  (B-14)  implies  that 


->  ( K,  0)  =  ^1-K  =  z/R 


so  that 


j  *1  =  r* 

J  ^  z  kj  (K,  0)  • 


(B-27) 


(B-28) 


Then  (B-22)  gives 


(B-29) 


|ki<K,z)E(K.z)H(K,z)  4n  Ik^K.O) 


Since  A  and  t  depend  only  on  K,  we  need  only  employ  the  magnitude  of  (B-16) 


in  the  sequel : 


|x-£|  =  K  E(K,  z)  . 


(B-30) 


Equation  (B-14),  (B-17),  (B-29)  and  (B-30)  are  equivalent  to  equations  (3-8) 
of  the  text. 


Finally,  we  relate  our  parameter  o  along  the  rays  to  the  corresponding 


arclength  parameter.  From  (B-7)  we  have 


dx  dx  dz  dz 

—  —  .  i  > 

—  •  —  + - =  K  +kj=n 

do  do  do  do 


(B-31) 


and  so 


'  A'A'. 


APPENDIX  C 


IHCHBOFF  DATA 


The  derivation  presented  in  Cohen  and  Bleistein  [1983]  applies  here  with 
the  constant  c  Green's  function  used  there  being  replaced  with  the  c(z) 
Green's  function  derived  in  Appendix  B.  Thus,  equation  (8)  of  that  paper 
may  be  recast  in  our  present  notation  as 


V* 


(z.  <o) 


R 


(x>  z;£/ w) 


(C-l) 


with  g  defined  by  equations  (3-d)  and  the  reflection  coefficient  R  being 
defined  by 


R  = 


y 

y  +  r. 


(C-2) 


where 


y  =  8  •  Vt  ,  Yj  =  sgn(Y> 
and  8  is  the  upward  normal  to  S.  Here  cx  is  defined  by  (31). 


(C-3) 


Since  g  has  the  form  given  in  (B-l) 


d  i 
Tn  * 


2im  8  •  Vt  g*  , 


( C-4) 


subject  to 


l?-£l  =  KE  (K.  z) 


( C-5) 


Hence  (B-l),  with  the  definition  of  y  given  in  (C-3)  and  that  of  S  given  by 
(29),  yields 


d  * 
TK  g 


2iwyS 


2iurc 

e 


(C-6) 


Thus  (26)  follows  from  the  form,  (25),  of  the  surface.  It  remains  to 
establish  the  detailed  calculation  of  y  given  in  (28). 


We  have: 


y  =  fi  •  Vt 

--OTT*  ' 


1  <h'.0,-l>  r„  91  ,  II  .  1 

-  ,(()  ■  l  °i  77  •  ®i  w  «  ®«  J 

11+(h') 


1  (h',0,-1) 

cToT 


\ 


l+(h ' ) 


-*  fh'G^-G^-G  ]  . 

—  L  K  3x  I  8z  z  J 


From  the  constraint  (C-5)  we  compute 


3K  _  x-£  dK 

dx  p(KEiK  '  dz 


KE 

ttit: 


(C-7) 


(C-8) 


and  then  (28)  follows  from  (C-7),  (C-8)  and  the  results  of  Appendix  A. 


STATION  ART  PHASE  CALCULATIONS 


Assuming  that  the  phase.  #(x).  has  a  single  simple  stationary  point  xg/ 


Vf(xs)  =  0 


det 


dx.dx. 
»  J 


(x  )  !«=  0 
-s 


(D-l) 


the  integral. 


I(X)  ~  J  eiX,(i>  A(x)dnx  ,  |x|  »  1  (D-2) 

can  be  evaluated  asymptotically  by  the  mnltidimensional  stationary  phase 
formula,  (see  Bleistein  [1984]  or  Bleistein  and  Handelsman  [1975]). 

I(X)  ~  f  -i* ” i  1  - - - exp(iXi  +  i  —  agn  X  sig  (#..)}  .  (D-3) 

lTrTJ  j|d.t«..Tf  *  ' 

Here  A(>  lg,  and  ««>■  denote  respectively  the  amplitude  A,  the  phase  i 
and  the  Hessian  matrix,  Of/dxjdxj)  evaluated  at  the  stationary  point, 
x  =  x$.  Further,  sig  (fjj)  denotes  the  signature  of  ( j ) >  i.e.  the  number 
of  positive  eigenvalues  minus  the  number  of  negative  eigenvalues. 

If  there  were  several  stationary  points,  (D-3)  would  be  replaced  by  a  sum 
over  the  contributions  from  these  several  points.  If  the  Hessian  matrix 
vanished,  (D-3)  would  have  to  be  replaced  by  a  more  general  result. 
However,  only  (D-3)  is  required  here. 

In  the  second  application  of  (D-5)  to  follow,  n  =  4.  Although  the 

evaluation  of  the  signature  of  a  symbolic  4x4  matrix  can  be  tedious  or 


impossible,  our  task  is  considerably  simplified  by  the  fact  that  our  Hessian 
matrices  have  the  special  form, 


i  = 


0  v  0 
P  o  a 
0  y  0 
u  0  6 


whence 


det  #  =  (ay  - 


v  )  ( p&  -  u  ) 


(D-4) 


(D-5) 


To  evaluate  the  signature,  we  need  to  determine  the  roots,  o,  of  the 
eigenvalue  equation. 


det  (*  -  ol)  =  0  ,  ( D-6 ) 

where  I  is  the  identity  matrix.  From  (D-5),  we  find  at  once: 

det(f  -  ol)  =  £  a*-  (a+y)  a  +  ay  -  v*  J  [  a*  -  (p+6)  a  +  p6  -  u*  J  .  (D-7) 

Thus, 

v*  >  ay  ,  u*  >  P6  =>  sig  (ijj)  =  0  •  (D-8) 


and  similarly, 

v  >  ay  ,  u  <  P&  =>  sig  (Ijj)  =  -2  ,  (D-9) 


etc. 

We  now  turn  to  the  specific  stationary  phase  calculations  in  the  text. 
We  shall  freely  use  the  results  of  Appendix  A  without  explicit  citation. 
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First  we  examine  the  phase  in  (12): 

f(£)  =  G(K,z)  -G(K'.z')  =  c(0)It(K.z)-t(I'.z')) 
where  G,  K,  and  K'  are  defined  by: 


G(K.z) 


(D-10) 


|x  -  =  IE  (K,  z)  *  EE  ,  (D-ll) 

|x'  -  sl  =  K'E(K'.Z')  *  K'E'  .  (D-12) 


Variables  x,  z.  x',  and  z'  £  z  are  viewed  as  fixed  for  this  calculation. 
Implicit  differentiation  of  (D-ll)  with  respect  to  £■  yields  (after 
appealing  to  Appendix  A), 


P 


(IE), 


a  k 


=  H 


a  k 


so 


a*  ti  -  «4  ac*  i4  -  Xj 

a?i  “Bai  '  a$.  "  K'E'H' 


Using  these  results  we  have 


1 


a*  ra*  as'l  5.  -  x.  f.  -  x: 

*_  n  1  1  1  1 

*i  a?.  K  3£.  "35.  E  ~  E' 

i  1  1. 


(D-13) 


To  obtain  onr  stationarity  condition  we  set  i ^  =  0  and  have 


h  -  *i  •; 


I  i  =  1,2  . 


Using  this  result  with  (D-ll)  and  (D-12)  leads  to 


K  =  K' 


(D-14) 


This  condition  says  that  the  critical  value  of  £  is  such  that  the  points 
r  *  (x,z)  and  r '  =  (x',z')  lie  on  the  same  ray  eaanating  from  (£.0).  See 
Fig.  1  above.  Ve  sybolically  denote  this  stationarity  choice  of  £  by  £$. 


We  now  proceed  to  calculate  the  Hessian  aatrix  (Ijj).  Differentiating 
expression  (D-13)  with  respect  to  £j  leads  to 


*ij  =  njr  =  8ij  &  ■  t\  ~  kikj  Kf  ■  eV] 


5  8ij e  -  kikj * 


fe  have  again  appealed  to  results  in  Appendix  A.  It  now  follows  that 


v.  *•-  * 


det  (f.j) 


e  -  ktr 


-kky 
1  2  1 


-kky 
1  2  * 


8  -  kjY 


a  a 

e  -  s  K  y 


=  e  (e  -  l  y) 


=  fL  _  i_l  Q_  _  J.  _  ** 

[£  E'J  [E  E’  Ira  E'B'JJ 


6  "id  Or -  if 
ft  “id  ft  “id 

E'  -  E  B*  -  B 
EE*  ISP 


[B-E  B'-E'l] 

LhT  “  E’H '  JJ 


Note  detdjj)  >  0  when  z  t  z'  since  E  and  H  are  z  integrals  with  positive 
integrands  (see  Appendix  A).  Similarly,  in  expanding  detlljj  -  «D  one 
finds  eigenvalues 


E'  -  E 


B*  -  B 


and  therefore  concludes  that  sig  (I,,)  =  2  sgn  (z'-z).  Note  that 


i  J  tigdij)  _  i  j  »gn  <z'-x>  .  A  ,gn  u._z)  . 


Using  these  resnlts  in  (D-3)  leads  to  the  following  asymptotic  (large  w) 
approximation  to  the  £  integral  in  (12): 


f.  V.  ■■ 

-  -A  »  -  u  ■ 


.  V- 

V-*  A'. 


* 

N 


v  .v 


c’(z)  k^K.O)  k  JK.z)  E(K,  z)  H ( K,  z ) 

in  c(0)  sgn  (z'-z)  (E  E'HH ' )  * _ 

w  [(E  -  E')  (B  -  B')]l/l 


We  now  turn  to  the  stationary  phase  evaluation  of  (32)  in  the  four 
variables  x,  y,  t,,  ^ .  Here  the  phase  is 

1  =  G(t.z)  -  G(K.z)  =  c(0)[t(I,z)  -  t(I.z)]  (D-15) 

subject  to  the  usual  constraints 

|x  -  =  K  E(K, z)  .  |x  -  £|  =  fE(l,z)  .  (D-16) 

This  stationary  phase  evaluation  is  soaewhat  aore  difficult  because  z 
depends  on  x: 

z  =  h(x)  ,  (D-17) 

Implicit  differentiation  in  x,  of  the  second  equation  in  (D-16)  leads  to 

=  (KE)r  —  +  IE  =g-^-  +  -5-h'  .  (D-l  8) 

p  dx  Z  dx  k( 


so 


dK  _ 

x-E 

Kh’ 

dz 

KEH 

k  H 
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(D-19) 


where  we  use  notations  like 


E  =  E(I,  z)  (D-20) 

for  quantities  which  depend  on  K  and  z.  The  remaining  partials  of  1  and  K 
are  simpler  since  z  does  not  depend  on  y,  or  y,  and  since  z  is  a 

completely  independent  variable: 


dK 

_  y~n 

3K 

5-x 

—  t 

dy 

KEH 

n 

KEH 

(D-21) 


These  results  allow  one  to  compute: 


i. 

X 


=  G  ii 
K  dx 


G 

z 


h '  (x) 


n*(z) 


h' 


p-y  _  n~y 

E  E 


(D-22) 


The  y,  r)  derivatives  give  rise  to  the  simple  stationarity  conditions: 


y  =  q  =  y  . 

Hence  the  constraints  (D-16)  reduce  to 

|x  -  5|  =  KE(K,  z)  .  |x  -  ?|  =  lE(K.z)  , 
which  can  be  rewritten  as 


(D-23 ) 


(D-24) 


x-$  =  pKE 


p  ■  sgn(x-!j) 


(D-25) 


P  •  Sgn(  x-1;) 


( D-26 ) 


x-$  =  pKE  , 


These  results  allow  us  to  simplify  the  remaining  partials  in  (D-22)  to 


i_  =  pi  +  k  h' 

x  j 


and 

l_  =  -  pi  +  pK 

For  stationarity,  we  must  have 


and  then  also 


p  =  p  =  -  sgn(h') 


K  =  K  =  k  h’  . 

i 

The  latter  equality  allows  us  to  determine  that 

i  =  I*:, l  ,  . 

^|l+h'  *  l|l+h'* 

Finally  (D-29)  allows  us  to  restate  (D-25),  (D-26)  as 

x  -  ?  =  -  I  E(I, z)  sgn(h') 

and 

x  -  5  =  -  KE(I,z)  sgn(h') 

At  the  stationary  point,  the  phase  is 


( D-27 ) 


(D-28) 


(D-29) 


(D-30) 


(D-31) 


(D-32) 


(D-33) 
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(D-34) 


|  =  G(K,  z)  -  G(K,z) 

s 


Farther,  a  number  of  terms  of  the  Hessian  matrix  vanish  becanse  of 
leaving  as  with 


ij 


*xx 

0 

0 

= 

0 

*yy 

0 

IH 

0 

*« 

0 

0 

V. 

0 

which  has  the  form  (D-4)  with 


-  ay  =  (#*5  -  !__  «K)|_  . 


and 


a2  -  pY  =  (|*_  -  I —  I  ) 

yn  yy  nn  s 


We  observe  that  (D-19)  and  the  stationarity  conditions  yield 


||  =  £  -  h'  lh'  1  =  JL  (l+h'1)  ,  p  =  -sgn(h') 

3x  s  H  H  H 


and  similarly 


3K  fi  dK 

*  ,"s  ’  h 


s  H(K,  z) 


(D-23 ) 


(D-3  5) 


(D-36) 


(D-37) 


(D-3  8) 


(D-39) 


These  facts  allow  to  evaluate  (D-36)  as 


Thus,  "near"  the  reflector,  z  =  z,  both  v1  >  ay,  u1  >  08  hold  and  from 
(D-8)  we  have  sig  (f^j)  =  0.  Finally,  applying  all  these  results  to  the 
integral  given  by  (32)  leads  to  the  approximation  (40). 


APPENDIX  E 


A  1-DIMENSIONAL  DISTRIBUTION 


In  (19)  we  were  lead  to  the  singular  form 


8  '  (t-t ' ) 
z-z ' 


z-z 


(E-l) 


This  requires  some  careful  interpretation  since  as  z  ->z',  z  ->  t '  and  hence 
the  form  is  very  singular  at  the  point  at  which  the  "action"  takes  place. 
We  here  first  show  that,  with 


p  =  |r-r ' |  sgn  (z-z ’)  , 


we  have 


lim  — 
z-z' 

z  ->  z  * 


n(z  * ) 

kj (K,  z ')  * 


(E-2) 


(E-3) 


Hm  ~T  =  lia  ~rr  =  c(z')  .  (E-4) 

T  *”X  T  “T 

z  ->  Z  1  T  ->T  ' 


To  establish  (E-3),  we  appeal  to  (B-16)  and  note  that  since  r  and  r* 


are  on  the  same  ray: 


r-r'  =  (x,  z)  -  (x\  z') 


Hence  as  z  -> z  * : 


Using  this  result  and  (E-2)  establishes  (E-3). 


To  obtain  (E-4)  we  use  (A-6)  to  get 


r-r  * 
z-z ' 


v  n  (z*) 

c<0)k|(E.s')  11 


z  ->  i  * 


Using  this  together  with  (E-3)  establishes  (E-4) 


Applying  (E-4)  we  have 


6(r-r ') 

4~  bix-x ') 
dr 


&[p/c(z')l  =  c(z')  6(p)  t 

c*(z')  -r—  6 < p)  =  c  (z’)6’(p) 
dp 


(E-6) 


This  result*  (E-6),  and  (E-3)  are  combined  in  the  main  text  to  further 
interpret  the  singular  expression  (E-l). 


Figure  1:  Two  ray«  (solid  lines)  connecting  an  arbitrary  surface  point  £ 
to  r  and  r'j  end  the  single  ray  eaianating  from  point  £s 


determined  by  stationary  phase. 


51  - 


uric  Lass  if  ied 


SECURITY  CLASSIFICATION  or  THIS  PAGE  (Whan  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 


1.  REPORT number 


CWP-021 


«.  TITLE  (and  Submit) 


READ  INSTRUCTIONS 
_ BEFC-R  COMPLETING  FORM 


S.  RECIPIENT'S  catalog  number 


S.  TYRE  or  REPORT  A  PERIOD  COVERED 


VELOCITY  INVERSION  USING  A  STRATIFIED  REFERENCE  Technical 


S.  PERFORMING  ORO.  REPORT  NUMBER 


7.  AUTMORfa; 


JACK  K.  COHEN  AND  FRANK  G.  HAG IN 


CONTRACT  OR  GRANT  NUmBERC«J 

N00014-84-K-0049 


t.  PERFORMING  ORGANIZATION  NAME  ANO  AOORESS 

Center  for  Wave  Phenomena 
Department  of  Mathematics 
Colorado  School  of  Mines,  Golden,  CO  80401 


II.  CONTROLLING  OFFICE  NAME  ANO  AOORESS 

Office  of  Naval  Research 
Arlington,  VA  22217 


I 


MONITORING  AGENCY  NAME  B  AOORESSfM  dlllata n<  ham  Centretllnd  OHIea) 


program  clement,  project,  task 

AREA  •  WORK  UNIT  NUMBERS 

NR  SRO-159/84APR20  (411) 


II.  REPORT  DATE 

January  25,  1985 


IS.  NUMBER  OF  PAGES 

51 


IS.  SECURITY  CLASS.  ( at  Mia  i 


IS.  DISTRIBUTION  STATEMENT  (»t  Mia  Report) 

This  document  has  been  approved  for  public  release  and  sale; 
its  distribution  is  ulimited. 


17.  DISTRIBUTION  STATEMENT  (ol  Ma  aAatracI  entered  In  Black  it,  II  dIHerent  ham  Report) 


20.  ABSTRACT  (Continue  an  taaataa  aide  II  nacaaaary  and  Idontlfr  kf  black  mmkaa) 


See  reverse  side 


roaa 
I  JAN  72 


EDITION  OF  I  NOV  SS  IS  OBSOLETE 
S/N  0102*0 14*  SSO I  I 


SECURITY  CLASSIFICATION  OF  THIS  PAOE  | 


i  Data  tnlerod) 


ABSTRACT 


The  purpose  of  this  work  is  to  present  an  inversion  algorithm  for 
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profile  in  realistic  earth  conditions.  The  basic  approach  follows  that  of 
the  original  Cohen  and  Bleistein  paper  [1979a]  in  that  high  frequency 
asymptotics  and  perturbation  methods  are  used.  However,  in  the  original 
paper  the  perturbation  was  relative  to  a  constant  reference  speed,  whereas 
the  current  work  uses  a  reference  speed  which  may  vary  with  depth.  This 
greatly  enhances  the  validity  of  the  perturbation  assumption  and  hence  the 
inversion  results.  On  the  other  hand,  the  new  algorithm  enjoys  the  same 
economies  and  stability  properties  of  the  original  algorithm,  making  it  very 
competitive  with  current  migration  schemes. 

Four  major  assumptions  are  made:  (i)  the  acoustic  wave  equation  is  an 
adequate  model,  (ii)  stacked  data  has  amplitude  information  worth  preserving 
fairly  accurately,  (iii)  the  actual  reflectivity  coefficients  can  be 
adequately  modeled  as  perturbations  from  a  continuous  reference  velocity 
which  depends  only  on  the  depth  variable,  and  (iv)  the  subsurface  can  be 
adequately  modeled  as  a  series  of  layers  with  jump  discontinuities  in  the 
velocity  (or  impedance)  at  these  layers. 

While  the  algorithm  is  particularly  suited  for  data  generated  by  a 
number  of  reflecting  surfaces,  its  validity  for  a  single  reflector  is 
demonstrated  by  applying  the  algorithm  to  Kirchhoff  data  for  a  quite  general 
surface . 

A  key  feature  of  the  approach  of  this  paper  is  the  repeated  application 
of  high  frequency  asymptotic  methods,*  both  in  obtaining  the  basic  integral 
equation  describing  the  unknown  velocity  correction,  and  in  the  inversion  of 
this  integral  equation.  Perhaps  a  noteworthy  feature  is  that  the  underlying 
integral  equation  is  in  the  form  of  a  generalized  Fourier  integral  equation,* 
and  the  method  for  its  (approximate)  inversion  may  prove  to  be  applicable  to 
a  wide  range  of  such  problems. 
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